#include "plugin.h"
#include <iostream>
#include <array>
#include <vector>
#include <unordered_map>
#include <memory>
#include <random>
#include <algorithm>
#include "tz_local_neuroseg.h"
#include "tz_trace_defs.h"
#include "tz_trace_utils.h"
#include "tz_locseg_chain_com.h"

namespace {
	typedef std::array<double, 3> Vec;
	Vec asVec(const Point& p) {
		Vec a;
		a[0]=p.x(); a[1]=p.y(); a[2]=p.z();
		return a;
	}
	Vec operator+(const Vec& a, const Vec& b) {
		Vec c;
		for(size_t i=0; i<b.size(); i++)
			c[i]=a[i]+b[i];
		return c;
	}
	Vec operator-(const Vec& a, const Vec& b) {
		Vec c;
		for(size_t i=0; i<b.size(); i++)
			c[i]=a[i]-b[i];
		return c;
	}
	Vec operator*(double a, const Vec& b) {
		Vec c;
		for(size_t i=0; i<b.size(); i++)
			c[i]=b[i]*a;
		return c;
	}
	Vec operator/(const Vec& b, double a) {
		Vec c;
		for(size_t i=0; i<b.size(); i++)
			c[i]=b[i]/a;
		return c;
	}
	double dot(const Vec& a, const Vec& b) {
		return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
	}
	double len(const Vec& a) {
		return sqrt(dot(a, a));
	}
	struct Seed {
		Edge edge;
		int index;
		Vec point;
		bool skipit;
	};
	struct Node {
		double time;
		double dist;
		int64_t index;
		size_t seed_i;
		bool skipit;
		bool closed;
	};

	//////////////////////

}


struct AlgFindNeutube : public AlgFind {
	double maxlen;
	double v0=0.0001;

	Edge prev_e;
	size_t prev_i;

	std::vector<Seed> seeds;
	std::vector<Node> store;
	std::vector<size_t> pq;
	std::unordered_map<int64_t, size_t> idx2node;

	double shell_dist_sum;
	double shell_dist_sum_delta;
	size_t nshell;

	//
	void setParams(const std::vector<ParamValue>& params) override {
		maxlen=params[0].getFloat();
	}
	std::vector<ParamValue> getParams() const override {
		std::vector<ParamValue> vl;
		vl.push_back(ParamValue{maxlen});
		return vl;
	}

	AlgFindNeutube() : AlgFind{}, maxlen{50.0} { }

	void config(Trace_Workspace*& tspace, Stack* stack, double* pos) {
		tspace=Locseg_Chain_Default_Trace_Workspace(nullptr, stack);
		if(!tspace->fit_workspace) {
			tspace->fit_workspace = New_Locseg_Fit_Workspace();
		}
		tspace->refit = true;
		tspace->tune_end = true;
		tspace->add_hit = true;
		tspace->min_score = 0.3;

		tspace->resolution[0]=cuberef<0>().xres;
		tspace->resolution[1]=cuberef<0>().yres;
		tspace->resolution[2]=cuberef<0>().zres;

		tspace->trace_range[0]=0;
		tspace->trace_range[3]=cuberef<0>().width-1;
		tspace->trace_range[1]=0;
		tspace->trace_range[4]=cuberef<0>().height-1;
		tspace->trace_range[2]=0;
		tspace->trace_range[5]=cuberef<0>().depth-1;
		if(pos[0]-maxlen/cuberef<0>().xres>tspace->trace_range[0])
			tspace->trace_range[0]=pos[0]-maxlen/cuberef<0>().xres;
		if(pos[0]+maxlen/cuberef<0>().xres<tspace->trace_range[3])
			tspace->trace_range[3]=pos[0]+maxlen/cuberef<0>().xres;
		if(pos[1]-maxlen/cuberef<0>().yres>tspace->trace_range[1])
			tspace->trace_range[1]=pos[1]-maxlen/cuberef<0>().yres;
		if(pos[1]+maxlen/cuberef<0>().yres<tspace->trace_range[4])
			tspace->trace_range[4]=pos[1]+maxlen/cuberef<0>().yres;
		if(pos[2]-maxlen/cuberef<0>().zres>tspace->trace_range[2])
			tspace->trace_range[2]=pos[2]-maxlen/cuberef<0>().zres;
		if(pos[2]+maxlen/cuberef<0>().zres<tspace->trace_range[5])
			tspace->trace_range[5]=pos[2]+maxlen/cuberef<0>().zres;
	}

	void destroy(Trace_Workspace*& tspace) {
		if (tspace) {
			if (tspace->fit_workspace) {
				Locseg_Fit_Workspace *fw =
					(Locseg_Fit_Workspace*) tspace->fit_workspace;
				fw->sws->mask = NULL;
				Kill_Locseg_Fit_Workspace(fw);
				tspace->fit_workspace = NULL;
			}
			Kill_Trace_Workspace(tspace);
			tspace = NULL;
		}
	}

	void vec2neutube(const Vec& v, double* p) {
		auto x=v[0]-cuberef<0>().origin[0];
		auto y=v[1]-cuberef<0>().origin[1];
		auto z=v[2]-cuberef<0>().origin[2];
		p[0]=x*cuberef<0>().rdir[0]+y*cuberef<0>().rdir[3]+z*cuberef<0>().rdir[6]-.5;
		p[1]=x*cuberef<0>().rdir[1]+y*cuberef<0>().rdir[4]+z*cuberef<0>().rdir[7]-.5;
		p[2]=x*cuberef<0>().rdir[2]+y*cuberef<0>().rdir[5]+z*cuberef<0>().rdir[8]-.5;
	}
	Vec neutube2vec(const double c[]) {
		auto x=c[0]+.5;
		auto y=c[1]+.5;
		auto z=c[2]+.5;
		auto px=x*cuberef<0>().dir[0]+y*cuberef<0>().dir[3]+z*cuberef<0>().dir[6]+cuberef<0>().origin[0];
		auto py=x*cuberef<0>().dir[1]+y*cuberef<0>().dir[4]+z*cuberef<0>().dir[7]+cuberef<0>().origin[1];
		auto pz=x*cuberef<0>().dir[2]+y*cuberef<0>().dir[5]+z*cuberef<0>().dir[8]+cuberef<0>().origin[2];
		return Vec{px, py, pz};
	}
	Point circle2point(const Geo3d_Circle& c) {
		auto x=c.center[0]+.5;
		auto y=c.center[1]+.5;
		auto z=c.center[2]+.5;
		auto px=x*cuberef<0>().dir[0]+y*cuberef<0>().dir[3]+z*cuberef<0>().dir[6]+cuberef<0>().origin[0];
		auto py=x*cuberef<0>().dir[1]+y*cuberef<0>().dir[4]+z*cuberef<0>().dir[7]+cuberef<0>().origin[1];
		auto pz=x*cuberef<0>().dir[2]+y*cuberef<0>().dir[5]+z*cuberef<0>().dir[8]+cuberef<0>().origin[2];
		Point p{px, py, pz, c.radius*(cuberef<0>().xres+cuberef<0>().yres)/2};
		return p;
	}
	Vec circle2vec(const Geo3d_Circle& c) {
		auto x=c.center[0]+.5;
		auto y=c.center[1]+.5;
		auto z=c.center[2]+.5;
		auto px=x*cuberef<0>().dir[0]+y*cuberef<0>().dir[3]+z*cuberef<0>().dir[6]+cuberef<0>().origin[0];
		auto py=x*cuberef<0>().dir[1]+y*cuberef<0>().dir[4]+z*cuberef<0>().dir[7]+cuberef<0>().origin[1];
		auto pz=x*cuberef<0>().dir[2]+y*cuberef<0>().dir[5]+z*cuberef<0>().dir[8]+cuberef<0>().origin[2];
		return Vec{px, py, pz};
	}

	void prepareCoreNodes(Edge e, int left, int right, int innerleft, int innerright) {
		auto nodeComp=[this](size_t a, size_t b)->bool {
			return store[a].time>store[b].time;
		};

		seeds.clear();
		store.clear();
		pq.clear();
		idx2node.clear();
		shell_dist_sum=0;
		shell_dist_sum_delta=0;
		nshell=0;

		double minres=cuberef<0>().xres;
		if(minres>cuberef<0>().yres)
			minres=cuberef<0>().yres;
		if(minres>cuberef<0>().zres)
			minres=cuberef<0>().zres;
		for(int k=left; k<right; k++) {
			auto pl=asVec(e.points()[k]);
			auto pr=asVec(e.points()[k+1]);
			auto dp=pr-pl;
			auto ldp=len(dp);
			auto nseg=lrint(3*ldp/minres);

			for(long i=0; i<=nseg; i++) {
				auto p=pl+i*dp/nseg;
				auto idx=xyzf2index<0>(p[0], p[1], p[2]);
				if(idx>=0) {
					index2xyzf<0>(idx, &p[0], &p[1], &p[2]);
					auto dp1=p-pl;
					auto r=dot(dp1, dp)/ldp;
					auto dist=sqrt(dot(dp1, dp1)-r*r);
					auto value=cube<0>(idx);
					auto time=dist/(value+v0);
					auto inode=idx2node.find(idx);
					Node* node=nullptr;
					if(inode!=idx2node.end()) {
						auto node=&store[inode->second];
						if(node->time<=time)
							node=nullptr;
					} else {
						auto ni=store.size();
						store.emplace_back();
						idx2node[idx]=ni;
						pq.push_back(ni);
						node=&store.back();
						node->seed_i=seeds.size();
						seeds.emplace_back();
					}
					if(node) {
						node->time=time;
						node->dist=dist;
						node->index=idx;
						//node->_i=//
						node->skipit=false;
						node->closed=false;
						seeds[node->seed_i].edge=e;
						seeds[node->seed_i].index=r*2>ldp?k+1:k;
						seeds[node->seed_i].point=r*dp/ldp+pl;
						seeds[node->seed_i].skipit=true;
					}
				}
			}
		}

		for(auto& n: store) {
			nshell++;
			shell_dist_sum+=n.dist;
			auto& s=seeds[n.seed_i];
			if(s.index>=innerleft && s.index<=innerright)
				s.skipit=false;
			else
				s.index=-1;
		}

		std::make_heap(pq.begin(), pq.end(), nodeComp);
	}

	Seed findSeed(Edge e, int cent, int left, int right, int innerleft, int innerright) {
		double maxres=cuberef<0>().xres;
		if(maxres<cuberef<0>().yres)
			maxres=cuberef<0>().yres;
		if(maxres<cuberef<0>().zres)
			maxres=cuberef<0>().zres;
		auto nodeComp=[this](size_t a, size_t b)->bool {
			return store[a].time>store[b].time;
		};

		if(pq.empty()) {
			prepareCoreNodes(e, left, right, innerleft, innerright);
		}

		Vec p0=asVec(e.points()[cent]);
		struct { int32_t dx, dy, dz; double dist; } all_neighbors [] = {
			{-1, 0, 0, xres<0>()},
			{1, 0, 0, xres<0>()},
			{0, -1, 0, yres<0>()},
			{0, 1, 0, yres<0>()},
			{0, 0, -1, zres<0>()},
			{0, 0, 1, zres<0>()}
		};

		size_t next_size=nshell+1000;
		while(!pq.empty() && sync(100)) {
			if(nshell>next_size) {
				next_size+=1000;
				//while(sync(100))
					//usleep(1000);
			}

			auto currentp=pq.front();
			std::pop_heap(pq.begin(), pq.end(), nodeComp); pq.pop_back();
			if(store[currentp].skipit)
				continue;

			Vec pp;
			auto current=store[currentp].index;
			index2xyzf<0>(current, &pp[0], &pp[1], &pp[2]);
			double dist0=len(pp-p0);
			if(dist0>maxlen) {
				if(seeds[store[currentp].seed_i].index>=0) {
					shell_dist_sum_delta+=store[currentp].dist;
					nshell--;
				}
				std::cerr<<"skip maxlen\n";
				continue;
			}
			if(store[currentp].dist>18) {
				if(seeds[store[currentp].seed_i].index>=0) {
					shell_dist_sum_delta+=store[currentp].dist;
					nshell--;
				}
				std::cerr<<"skip dist\n";
				continue;
			}

			int32_t x, y, z;
			index2xyzi<0>(current, &x, &y, &z);
			std::cerr<<"N: "<<nshell<<" "<<(shell_dist_sum-shell_dist_sum_delta)/nshell<<" "<<store[currentp].dist<<" "<<store[currentp].time<<" "<<x<<" "<<y<<" "<<z<<"\n";

			if((shell_dist_sum-shell_dist_sum_delta)/nshell>maxres && store[currentp].dist>(shell_dist_sum-shell_dist_sum_delta)/nshell*2.0) {
				if(!seeds[store[currentp].seed_i].skipit) {
					Seed seed;
					seed.edge=seeds[store[currentp].seed_i].edge;
					seed.index=seeds[store[currentp].seed_i].index;
					index2xyzf<0>(current, &seed.point[0], &seed.point[1], &seed.point[2]);
					seed.skipit=false;

					shell_dist_sum_delta+=store[currentp].dist;
					nshell--;
					std::cerr<<"new seed\n";
					return seed;
				} else {
					std::cerr<<"skip seed "<<store[currentp].seed_i<<" "<<seeds[store[currentp].seed_i].index<<"\n";
				}
			}
			/*
				*/

			store[currentp].closed=true;
			if(seeds[store[currentp].seed_i].index>=0) {
				shell_dist_sum_delta+=store[currentp].dist;
				nshell--;
			}

			auto value=cube<0>(current);
			for(auto& delta: all_neighbors) {
				int32_t neighbor_x=x+delta.dx;
				int32_t neighbor_y=y+delta.dy;
				int32_t neighbor_z=z+delta.dz;
				if(neighbor_x<0 || neighbor_x>=width<0>())
					continue;
				if(neighbor_y<0 || neighbor_y>=height<0>())
					continue;
				if(neighbor_z<0 || neighbor_z>=depth<0>())
					continue;

				Vec np;
				int64_t neighbor=xyzi2index<0>(neighbor_x, neighbor_y, neighbor_z);
				index2xyzf<0>(neighbor, &np[0], &np[1], &np[2]);

				auto dist=len(np-seeds[store[currentp].seed_i].point);
				auto time=store[currentp].time+delta.dist/(value+v0);
				auto neighbori=idx2node.find(neighbor);
				if(neighbori!=idx2node.end()) {
					auto nn=neighbori->second;
					if(time<store[nn].time) {
						if(store[nn].closed) {
							store[nn].closed=false;
							store[nn].time=time;
							store[nn].dist=dist;
							store[nn].seed_i=store[currentp].seed_i;
							pq.push_back(nn); std::push_heap(pq.begin(), pq.end(), nodeComp);
							if(seeds[store[nn].seed_i].index>=0) {
								nshell++;
								shell_dist_sum+=dist;
							}
						} else {
							store[nn].skipit=true;
							if(seeds[store[nn].seed_i].index>=0) {
								nshell--;
								shell_dist_sum_delta+=store[nn].dist;
							}
							nn=store.size();
							store.emplace_back();
							store[nn].time=time;
							store[nn].dist=dist;
							store[nn].skipit=false;
							store[nn].closed=false;
							store[nn].index=neighbor;
							store[nn].seed_i=store[currentp].seed_i;
							idx2node[neighbor]=nn;
							pq.push_back(nn); std::push_heap(pq.begin(), pq.end(), nodeComp);
							if(seeds[store[nn].seed_i].index>=0) {
								nshell++;
								shell_dist_sum+=dist;
							}
						}
					}
				} else {
					auto nn=store.size();
					store.emplace_back();
					store[nn].time=time;
					store[nn].dist=dist;
					store[nn].skipit=false;
					store[nn].closed=false;
					store[nn].index=neighbor;
					store[nn].seed_i=store[currentp].seed_i;
					idx2node[neighbor]=nn;
					pq.push_back(nn); std::push_heap(pq.begin(), pq.end(), nodeComp);
					if(seeds[store[nn].seed_i].index>=0) {
						nshell++;
						shell_dist_sum+=dist;
					}
				}
			}
		}
		Seed seed{Edge{}};
		seed.skipit=true;
		return seed;
	}

	Seed findSeed(const Position& p0) {
		if(prev_e!=p0.edge || prev_i!=p0.index) {
			pq.clear();
			prev_e=p0.edge;
			prev_i=p0.index;
		}
		if(p0.edge) {
			size_t l, r;
			for(r=p0.index; r<p0.edge.points().size(); r++) {
				if(p0.edge.points()[r].distTo(p0.point)>maxlen)
					break;
			}
			r--;
			for(l=p0.index; l-->0; ) {
				if(p0.edge.points()[l].distTo(p0.point)>maxlen)
					break;
			}
			l++;
			if(p0.index==0) {
				return findSeed(p0.edge, p0.index, l, r, p0.index, p0.index);
			} else if(p0.index==p0.edge.points().size()-1) {
				return findSeed(p0.edge, p0.index, l, r, p0.index, p0.index);
			} else {
				return findSeed(p0.edge, p0.index, l, r, l, r);
			}
		} else {
			Seed seed;
			seed.edge=Edge{};
			seed.index=-1;
			seed.point=asVec(p0.point);
			seed.skipit=false;
			return seed;
		}
	}
	// ////// // // // // // // // 

	std::vector<size_t> integrateResultX(Edge edge, size_t index, const std::vector<Point>& points) {
		std::vector<size_t> idxes{};
		if(points.size()<2)
			return idxes;
		if(!edge) {
			for(size_t i=0; i<points.size(); i++)
				idxes.push_back(i);
			return idxes;
		}
		std::vector<size_t> idx_map(points.size(), SIZE_MAX);
		std::vector<size_t> idx_rmap(points.size(), SIZE_MAX);
		int mi;
		double pdist=INFINITY;
		for(size_t i=0; i<points.size(); i++) {
			auto d=points[i].distTo(edge.points()[index]);
			if(d<pdist) {
				pdist=d;
				mi=i;
			}
		}
		idx_map[mi]=index;
		idx_rmap[mi]=index;
		std::cerr<<"PDIST: "<<pdist<<" "<<"\n";
		double mdist00, mdist01, mdist10, mdist11;
		{
			size_t j0=index;
			double sumdist=0;
			for(size_t i=mi+1; i<points.size(); i++) {
				double mdist=INFINITY;
				double j1;
				for(size_t j=j0; j<j0+5 && j<edge.points().size(); j++) {
					auto d=points[i].distTo(edge.points()[j]);
					if(d<mdist) {
						mdist=d;
						j1=j;
					}
				}
				sumdist+=mdist;
				idx_map[i]=j1;
				j0=j1;
			}
			mdist11=(sumdist+pdist)/(points.size()-mi);
			std::cerr<<"MDIST11: "<<mdist11<<" "<<points.size()-mi<<"\n";
		}
		{
			size_t j0=index;
			double sumdist=0;
			for(size_t i=mi; i-->0; ) {
				double mdist=INFINITY;
				double j1;
				for(size_t j=j0+1; j+4>j0 && j-->0; ) {
					auto d=points[i].distTo(edge.points()[j]);
					if(d<mdist) {
						mdist=d;
						j1=j;
					}
				}
				sumdist+=mdist;
				idx_map[i]=j1;
				j0=j1;
			}
			mdist00=(sumdist+pdist)/(mi+1);
			std::cerr<<"MDIST00: "<<mdist00<<" "<<mi+1<<"\n";
		}
		{
			size_t j0=index;
			double sumdist=0;
			for(size_t i=mi+1; i<points.size(); i++) {
				double mdist=INFINITY;
				double j1;
				for(size_t j=j0+1; j+4>j0 && j-->0; ) {
					auto d=points[i].distTo(edge.points()[j]);
					if(d<mdist) {
						mdist=d;
						j1=j;
					}
				}
				sumdist+=mdist;
				idx_rmap[i]=j1;
				j0=j1;
			}
			mdist10=(sumdist+pdist)/(points.size()-mi);
			std::cerr<<"MDIST10: "<<mdist10<<" "<<points.size()-mi<<"\n";
		}
		{
			size_t j0=index;
			double sumdist=0;
			for(size_t i=mi; i-->0; ) {
				double mdist=INFINITY;
				double j1;
				for(size_t j=j0; j<j0+5 && j<edge.points().size(); j++) {
					auto d=points[i].distTo(edge.points()[j]);
					if(d<mdist) {
						mdist=d;
						j1=j;
					}
				}
				sumdist+=mdist;
				idx_rmap[i]=j1;
				j0=j1;
			}
			mdist01=(sumdist+pdist)/(mi+1);
			std::cerr<<"MDIST01: "<<mdist01<<" "<<mi+1<<"\n";
		}

		bool left;
		if(index==0) {
			if(mdist01>mdist11) {
				left=true;
			} else {
				left=false;
			}
		} else if(index==edge.points().size()-1) {
			if(mdist10<mdist00) {
				left=true;
			} else {
				left=false;
			}
		} else {
			if(mdist00*mdist01<mdist10*mdist11) {
				left=false;
			} else {
				left=true;
			}
		}
		bool usem=true;
		if(pdist<edge.points()[index].r())
			usem=false;

		if(!left) {
			if(usem)
				idxes.push_back(mi);
			for(size_t i=mi+1; i<points.size(); i++)
				idxes.push_back(i);
		} else {
			if(usem)
				idxes.push_back(mi);
			for(int i=mi; i-->0; )
				idxes.push_back(i);
		}

	}
	std::vector<size_t> integrateResult(Edge edge, int index, const std::vector<Point>& points) {
		std::vector<size_t> idxes{};
		if(points.size()<2)
			return idxes;
		if(!edge) {
			for(size_t i=0; i<points.size(); i++)
				idxes.push_back(i);
			return idxes;
		}

		typedef std::tuple<size_t, size_t, double, bool> Seg;
		std::vector<Seg> segments;
		bool inside=false;
		for(size_t i=0; i<points.size(); i++) {
			auto v=asVec(points[i]);
			auto idx=xyzf2index<0>(v[0], v[1], v[2]);
			bool hit=false;
			if(idx>=0) {
				auto nodei=idx2node.find(idx);
				if(nodei!=idx2node.end()) {
					if(store[nodei->second].closed)
						hit=true;
				}
			}
			if(!hit) {
				if(!inside) {
					inside=true;
					segments.emplace_back();
					std::get<0>(segments.back())=i;
				}
				std::get<1>(segments.back())=i;
			} else {
				inside=false;
			}
		}
		for(auto& seg: segments) {
			auto dl=points[std::get<0>(seg)].distTo(edge.points()[index]);
			auto dr=points[std::get<1>(seg)].distTo(edge.points()[index]);
			if(dl<dr) {
				std::get<2>(seg)=dl;
				std::get<3>(seg)=false;
			} else {
				std::get<2>(seg)=dr;
				std::get<3>(seg)=true;
			}
		}
		std::sort(segments.begin(), segments.end(), [](const Seg& a, const Seg& b) { return std::get<2>(a)<std::get<2>(b); });
		if(segments.size()<=0)
			return idxes;
		auto& seg=segments[0];
		auto maxl=std::get<0>(seg);
		auto maxr=std::get<1>(seg);
		if(maxr-maxl+1<2)
			return idxes;

		//
		//
		//
		//, maxlen=0, maxl;
			//std::cerr<<"hit "<<hit<<"\n";
		auto pc=asVec(edge.points()[index]);
		if(!std::get<3>(seg)) {
			while(maxl>0) {
				auto p1=asVec(points[maxl-1]);
				auto p0=asVec(points[maxl]);
				auto dp0c=pc-p0;
				auto dp01=p1-p0;
				auto r=dot(dp01, dp0c)/dot(dp01, dp01);
				if(r<=1)
					break;
				maxl--;
			}
			for(size_t i=maxl; i<=maxr; i++)
				idxes.push_back(i);
		} else {
			while(maxr+1<points.size()) {
				auto p1=asVec(points[maxr+1]);
				auto p0=asVec(points[maxr]);
				auto dp0c=pc-p0;
				auto dp01=p1-p0;
				auto r=dot(dp01, dp0c)/dot(dp01, dp01);
				if(r<=1)
					break;
				maxr++;
			}
			for(size_t i=maxr+1; i-->maxl; )
				idxes.push_back(i);
		}
		return idxes;
	}


	void find(const Position& p0, Graph gr, Path* pth) override {
		pth->edge1=Edge{};
		pth->index1=-1;
		pth->points.clear();

		Stack stack;
		stack.kind=cuberef<0>().is16bit?GREY16:GREY;
		stack.width=cuberef<0>().widthAdj;
		stack.height=cuberef<0>().heightAdj;
		stack.depth=cuberef<0>().depthAdj;
		stack.text=nullptr;
		stack.array=const_cast<uint8*>(reinterpret_cast<const uint8*>(cuberef<0>().data));

		/*
			while(sync(100))
			usleep(1000);
			*/


		do {
			Trace_Workspace* tspace=nullptr;

			auto seed=findSeed(p0);
			if(seed.skipit) {
				printMessage("Null seed\n");
				return;
			}

			printMessage("seed: ", seed.point[0],' ', seed.point[1],' ', seed.point[2]);

			double pos[3];
			vec2neutube(seed.point, pos);
			config(tspace, &stack, pos);

			Local_Neuroseg *locseg = New_Local_Neuroseg();
			Locseg_Fit_Workspace *ws = (Locseg_Fit_Workspace*) tspace->fit_workspace;

			Set_Neuroseg(&(locseg->seg), 3.0, 0.0, 11.0, TZ_PI_4, 0.0, 0.0, 0.0, 1.0);
			Set_Neuroseg_Position(locseg, pos, NEUROSEG_CENTER);
			Local_Neuroseg_Optimize_W(locseg, &stack, 1.0, 1, ws);
			double optpos[3];
			Get_Neuroseg_Position(locseg,  optpos);

			Trace_Record *tr = New_Trace_Record();
			tr->mask = ZERO_BIT_MASK;
			Trace_Record_Set_Fix_Point(tr, 0.0);
			Trace_Record_Set_Direction(tr, DL_FORWARD);
			Locseg_Node *p = Make_Locseg_Node(locseg, tr);
			Locseg_Chain *locseg_chain = Make_Locseg_Chain(p);

			Trace_Workspace_Set_Trace_Status(tspace, TRACE_NORMAL, TRACE_NORMAL);
			Trace_Locseg(&stack, 1.0, locseg_chain, tspace);
			Locseg_Chain_Remove_Overlap_Ends(locseg_chain);
			Locseg_Chain_Remove_Turn_Ends(locseg_chain, 1.0);

			int n = 0;
			Geo3d_Circle *circles=Locseg_Chain_To_Geo3d_Circle_Array(locseg_chain, NULL, &n);
			Kill_Locseg_Chain(locseg_chain);

			std::vector<Point> points{};
			bool inside=false;
			int inside_i;
			int f=0, t=0, ftl=0;
			for(int i = 0; i < n; ++i) {
				auto p=circle2point(circles[i]);
				auto d=p.distTo(p0.point);
				if(d<maxlen) {
					if(!inside) {
						inside=true;
						inside_i=i;
					}
					if(i-inside_i+1>ftl) {
						f=inside_i;
						t=i;
						ftl=t-f+1;
					}
				} else {
					inside=false;
				}
			}
			for(int i=f; i<=t; i++)
				points.push_back(circle2point(circles[i]));
			//auto seed_cur=seed_i;
			{
				auto nodeComp=[this](size_t a, size_t b)->bool {
					return store[a].time>store[b].time;
				};
				//auto p=neutube2vec(optpos);
				//auto idx=xyzf2index<0>(p[0], p[1], p[2]);
				auto idx=xyzf2index<0>(seed.point[0], seed.point[1], seed.point[2]);
				if(idx>=0) {
					auto inode=idx2node.find(idx);
					Node* node=nullptr;
					if(inode!=idx2node.end()) {
						node=&store[inode->second];
						if(node->time<=0.0)
							throwError("not possible");
						if(node->closed)
							throwError("not possible");
						node->time=0.0;
						node->dist=0.0;
						node->index=idx;
						node->seed_i=seeds.size();
						seeds.emplace_back();
						seeds.back().edge=Edge{};
						seeds.back().index=-1;
						seeds.back().point=seed.point;
						seeds.back().skipit=true;
						node->skipit=false;
						node->closed=false;
						pq.push_back(inode->second); std::push_heap(pq.begin(), pq.end(), nodeComp);
						//if(-1>=0) {
							//nshell++;
							//shell_dist_sum+=node->dist;
						//}
						std::cerr<<"insert new core\n";
					}
				}
			}
			/*
				for(int j=0; j<n; ++j) {
				auto p=circle2vec(circles[j]);
				auto idx=xyzf2index<0>(p[0], p[1], p[2]);
				for(size_t i=seed_i+1; i<seeds.size(); i++) {
				if(len(seeds[i].point-p)<circles[j].radius)
				seeds[i].skipit=true;
				}
				}
				seed_i++;
				while(seed_i<seeds.size() && seeds[seed_i].skipit)
				seed_i++;
				*/
			if (circles != NULL)
				free(circles);
			destroy(tspace);
			auto idxes=integrateResult(seed.edge, seed.index, points);

			if(idxes.size()>=2) {
				pth->edge0=seed.edge;
				pth->index0=seed.index;
				if(pth->edge0) {
					pth->points.push_back(pth->edge0.points()[pth->index0]);
				}
				//
				for(auto i: idxes) {
					pth->points.push_back(points[i]);
				}
				break;
			}
			//seed_i++;
		} while(true);
	}

	std::vector<int32_t> getUpdate() const override {
		std::vector<int32_t> data{};
		/*
			for(size_t i=seed_i; i<seeds.size(); i++) {
			auto& seed=seeds[i];
			Point p{seed.point[0], seed.point[1], seed.point[2]};
			data.push_back(p._x);
			data.push_back(p._y);
			data.push_back(p._z);
			}
			std::cerr<<"seeds size: "<<data.size()<<"\n";
			*/
		for(size_t i=0; i<pq.size(); i++) {
			auto& n=store[pq[i]];
			if(!n.skipit) {
				double xf, yf, zf;
				index2xyzf<0>(n.index, &xf, &yf, &zf);
				Point p{xf, yf, zf};
				data.push_back(p._x);
				data.push_back(p._y);
				data.push_back(p._z);
			}
		}
		return data;
	}
	void positionChanged() override {
		store.clear();
		seeds.clear();
		pq.clear();
		idx2node.clear();

		shell_dist_sum=0;
		shell_dist_sum_delta=0;
		nshell=0;
	}

	static AlgFind* creater() {
		return new AlgFindNeutube{};
	}
};

static void neutube_load(std::vector<AlgorithmInfo*>& algs, std::vector<GraphReaderInfo*>& readers) {
	auto infofind=new AlgInfo<AlgFind>{"NeuTube", 80,
		"Find neurites with neutube.",
		{
			ParamSpec{"Max length", "Maximum length.", 1.0, 255.0, 50.0, "µm"},
		},
		AlgFindNeutube::creater};
	algs.push_back(infofind);
}

static void neutube_unload() {
}

FNT_PLUGIN("NeuTube", &neutube_load, &neutube_unload, "Find neurites with NeuTube.");

